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Employing large-scale quantum Monte Carlo simulations, we reveal the full phase diagram of the 
extended Hubbard model of hard-core bosons on the pyrochlore lattice with partial fillings. When 
the inter-site repulsion is dominant, the system is in a cluster Mott insulator phase with integer 
number of bosons localized inside the tetrahedral nnits of the pyrochlore lattice. We show that the 
full phase diagram contains three cluster Mott insulator phases with 1/4, 1/2 and 3/4 boson fillings, 
respectively. We further demonstrate that all three cluster Mott insulators are Coulomb liquid phases 
and its low-energy property is described by the emergent compact U{1) quantum electrodynamics 
(QED). Besides measuring the specific heat and entropy of the cluster Mott insulators, we investigate 
the correlation function of the emergent electric field and verify it is consistent with the compact 
U{1) quantum electrodyamics description. Our results shed lights on magnetic properties of various 
pyrochlore systems as well as the charge physics of the cluster magnets. 


Introduction. Cluster Mott insulators (CMIs) and re¬ 
lated materials are new class of physical systems that 
provides a new ground to look for emergent and exotic 
phenomena m- Unlike in the conventional Mott insu¬ 
lators, the electrons in CMIs are localized inside certain 
cluster units rather than on the lattice sites. It is the 
inter-site electron interaction that leads to the uncon¬ 
ventional electron localization in CMIs and causes the 
electron motion inside neighboring clusters to be corre¬ 
lated. This correlated electron motion allows the elec¬ 
tron charges to fluctuate quantum mechanically in the 
positional space at energies below the Mott gap. Such 
a sub-gap charge fluctuation is one of the key features 
of CMIs and is responsible for various emergent phenom¬ 
ena [T] . For instance, in the electron CMI on a pyrochlore 
lattice where the electrons are localized on the tetrahe¬ 
dral clusters, it was shown that [1], due to the sub-gap 
charge fluctuation, the charge sector is in a Coulomb liq¬ 
uid phase with emergent gapless 17(1) gauge photon and 
fractional charge excitation. 

The existing proposals for CMIs focus on the electron 
systems as the real materials are electronic systems. Due 
to the fermionic statistics, however, the relevant theoret¬ 
ical models often suffer from a fermion sign problem and 
are prohibited from controlled investigation via quan¬ 
tum Monte Carlo (QMC) simulations. Furthermore, the 
electron models are complicated by the presence of the 
spin carrying excitations. Nevertheless, the hallmark of 
the CMIs - particles being localized on the cluster units 
with sub-gap fluctuations - is independent of the particle 
statistics. Therefore, in this paper, we turn our attention 
to bosons and study the hard-core boson models on a 
pyrochlore lattice with partial fillings. We show that the 



FIG. 1. (Color online), (a) Snapshot of the boson config¬ 
uration on the pyrochlore lattice at CMI-1/4, viewed on a 
Kagome plane. Removing one boson creates two defect tetra- 
hedra (marked with dashed circles). The intensity of the con¬ 
tour background stands for the l/r'* decay of the electric field 
correlation for the Coulomb liquid phase, where the electric 
field component normal to the Kagome plane is chosen. The 
vectors ai, 2,3 are the primitive vectors of the FCC Bravais lat¬ 
tice. (b) The sub-gap boson density fluctuations of the CMI- 
1/4, happening as the collective tunneling of three bosons 
around the perimeter of the elementary hexagon. 


bosonic pyrochlore system supports the Coulomb liquid 
phase with emergent 17(1) gauge structure in the bosonic 
CMI. The emergent and exotic physics in the charge sec¬ 
tor of the electronic CMI is fully retained in the corre¬ 
sponding bosonic CMI on the same lattice. 
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Model. Our model is defined on the pyrochlore lattice 
with the Hamiltonian 

i? = ^ [-t{blbj + h.c.) + Vnitij] - (1) 

(*J> i 

where b^ (b) is the hard-core boson creation (annihila¬ 
tion) operator. By definition, we exclude the double oc¬ 
cupation of the bosons on a single lattice site, t {V) is 
the nearest neighbor boson hopping (repulsion), and we 
set y = 1 as the energy unit throughout the paper. 

We first clarify the nomenclatures associated with this 
model and its ground state in the literature [ 11 ]. The 
model in Eq. was originally proposed as the spin-1/2 
XXZ model on a pyrochlore lattice in Ref. H if the hard¬ 
core boson variables are transformed into spin operators 
via b\ = Sf+iSj and rii —1/2 = Sf. The hard-core boson 
model at 1/2 boson filling gives a Coulomb liquid CMI 
phase when the inter-site repulsion V is strong enough to 
localize the bosons in the tetrahedra with two bosons per 
tetrahedron HD. This insulating Coulomb liquid phase 
at 1/2 filling is referred as “CMI-1/2” (see Fig. in 
this letter. In the spin language, the “CMI-1/2” phase is 
the well-known quantum spin ice phase whose low-energy 
property is described by the compact U{1) gauge field 
coupled with gapped bosonic spinons that carry quan¬ 
tum nubmer = 1/2 (THUj, and the boson localization 
condition with two bosons in each tetrahedron is essen¬ 
tially the “two-up two-down” spin ice rule |4]. 

In this letter, we carry out a large-scale world-line 
QMC simulation to investigate the full phase diagram of 
Eq.[T] in a more general parameter space. We vary both 
t/V and p/V to access different interaction strength as 
well as the boson fillings other than the 1/2 filling. We 
reveal two more CMIs at 1/4 and 3/4 fillings (see the 
schematic plot of the CMI at 1/4 filling in Fig. [^(a) and 
the phase diagram in Fig. and verify both of them are 
in Coulomb liquid phase with emergent compact QED 
description of the low-energy property. 

Numerical methods. We employ large-scale worm- 
type QMC technique to simulate the Hamiltonian in 
Eq.[T] on the pyrochlore lattice. The simulations are 
performed in grand canonical ensemble. The worm al¬ 
gorithm is based on the path-integral formulation of 
the quantum many-body action and works in continu¬ 
ous imaginary time without the Trotter error [SHE]- 
In the simulation, we choose the linear lattice size L in 
the range of 4 < L < 16, larger than those in previous 
studies mi. The inverse temperatures goes as large as 
/3 = lOOL, while simulations at even lower temperatures 
were also performed for some cases. 

To determine the phase boundary between the super- 
fluid and CMIs, we measure both the superfluid density 
Ps = + ITaa + s-iid the boson parti¬ 

cle density p = (iV)/(4L^), where {N) is the expectation 
value of the number of bosons and Wg ,., with i = 1,2, 3, 



FIG. 2. (Color online). Phase diagram of the hard-core boson 
model on a pyrochlore lattice. CMIs are at 1/4, 1/2 and 3/4 
boson fillings, respectively. Phase boundaries are determined 
in QMC simulations with system sizes 4 x L®. The physical 
properties along the (blue) dashed paths and at the (magenta) 
solid dot are plotted in Figs. [^ [^ and 

are the winding numbers [16] along the three primitive 
vector directions shown in Fig. &■ 

Phase diagram. The full phase diagram obtained 
from our simulations is presented in Fig. When t/V 
is small, two CMIs at 1/4 and 3/4 boson fillings in addi¬ 
tion to the known CMI at 1/2 filling are found. Here, the 
“CMI-1/4” (“CMI-3/4”) contains one (three) bosons in 
each tetrahedral cluster. When t/V is large, superfluid 
phase is obtained. 

In Figs. and b, we depict the boson density p and 
the superfluid density ps as a function of interaction t/V 
for a fixed chemical potential p/V = 1, i.e., the system 
is inside the CMI-1/4 when t/V is small. One sees in 
the superfiuid, ps is finite while p is strongly fluctuating. 
The phase boundary of the superfluid is signalled by the 
vanishing of ps. In Figs. H and P, we depict the boson 
density p and the superfluid density ps as a function of 
chemical potential p/V for a fixed interaction strength 
t/V = 0.04, i.e., the system is inside the CMI-1/4 and 
CMI-1/2 when p/V is within the Mott regime. In the 
Mott regime, the boson density is fixed while the super¬ 
fiuid density is zero. A fixed boson density implies that 
the compressibility dp/dp is zero and the system is hence 
Mott insulating. The Mott insulating phase is terminated 
at the places where the compressibility is non-zero. In 
all the cases that we simulated in Fig. and Fig. the 
phase boundary of the Mott insulators coincides with the 
phase boundary of the superfiuid. Therefore, there is no 
other intermediate phase between the superfiuid and the 
Mott insulators. 

One can understand the phase diagram in Fig. [^start¬ 
ing from the limit t —>■ 0. The superfiuid extends down to 
the t = 0 at even integer values of p/V, where there is no 
energy cost to the addition or removal of one bosons onto 
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FIG. 3. (Color online), (a) Boson density p and (b) superflnid 
density ps as a function of interaction tjV along the dashed 
horizontal path in Fig. with pjV = 1. (c) Boson density 
p and (d) superflnid density ps as a function of chemical po¬ 
tential p,/V along the dashed vertical path in Fig. with 
t/V = 0.04. 


the tetrahedral unit, and hence superfluidity occurs at ar¬ 
bitrarily small hopping t. Notice that adding (removing) 
one boson to the system creates two particle-like (hole¬ 
like) defect tetrahedra, as shown in Fig. & for the CMI- 
1/4. Here the particle-like (hole-like) defect tetrahedron 
refers to the tetrahedron that contains one more (less) 
boson than the boson number demanded by the Mott 
localization. The Mott lobe at the finite hopping t is un¬ 
derstood as the condensation of the mobile particle-like 
or hole-like defect tetrahedra. As t is increased from zero 
for a fixed chemical potential inside the Mott regime, the 
kinetic energy gain of the particle-like or hole-like defect 
tetrahedra increases. When the kinetic energy gain even¬ 
tually balances the interaction energy cost, the particle¬ 
like or hole-like excitations would like to condense, giv¬ 
ing rise to superfluidity. Increasing t reduces the chemi¬ 
cal potential width of the Mott regime and leads to the 
lobe-like phase boundary in Fig. 

At first looking, our phase diagram and the Mott lobes 
are quite similar to the ones in the conventional boson 
Hubbard model but there are key distinctions. In 
conventional Mott insulators, the boson is at integer fill¬ 
ings and the on-site Hubbard interaction freezes the bo¬ 
son motion and localizes the bosons on the lattice sites. 
In contrast, the CMIs in Fig. are at fractional fill¬ 
ings and it is the inter-site repulsion that localizes them, 
moreover, the bosons are localized inside the tetrahedral 
clusters of the pyrochlore lattice. For the CMI-1/4 (CMI- 
3/4) in our phase diagram, one boson (three bosons) is 
(are) localized on each tetrahedral unit. 

More important difference between the conventional 
Mott insulators and the CMIs in Fig. lies in the low- 
energy density fluctuations. In conventional Mott insula¬ 


tors, the low-lying density fluctuation that conserves the 
total boson number is the particle-hole excitation and 
costs an energy of the order of the Mott gap. In the CMIs, 
in addition to the similar type of particle-hole excitation 
as in the conventional Mott insulators, there exist boson 
density fluctuations far below the Mott gap. As shown 
in Fig. [^, for the CMI-1/4, this sub-gap boson density 
fluctuation arises from the collective tunnelling of three 
bosons around the perimeter of the elementary hexagon 
on the pyrochlore lattice. This collective tunnelling is 
a third order process of the nearest neighbor hopping 
t in the Mott regime and costs an energy of 0{t^/V'^). 
This process preserves the center of mass of the three 
bosons and is hence below the Mott gap. Like in the 
CMI-1/2 (or quantum spin ice in the spin language) [1], 
it is the sub-gap density fluctuation that gives rise to the 
Coulomb liquid phase with emergent compact 17(1) QED 
description at low energies for both the CMI-1/4 and the 
CMI-3/4 [HUH]. 

Coulomb liquid phase of the CMI-1/4. Now we focus 
on the CMI-1/4 and provide two numerical evidences, i.e. 
the correlation function of the emergent electric field, as 
well as the temperature dependence of the heat capacity 
and entropy, to support the compact 17(1) QED nature 
of the Coulomb liquid phase. 

In the compact QED description of the Coulomb liq¬ 
uid phase mini, the emergent electric field operator 
is related to the boson density operator as E(r) = 
~ l/4)ei, where the sum is over the four lattice 
sites on the tetrahedron centered at r and is the unit 
vector that points inwards (outwards) from r if r belongs 
to the up-pointing (down-pointing) tetrahedron. Using 
the compact QED action appropriate for the low-energy 
properties of the Coulomb liquid phase |3], it is readily 
to evaluate the equal-time electric field correlation and 
obtain 


{E^{k,T)E,(k',T)) 


- ki^k^)S(k^ 


■k') 


{2TTflKc flVp\k\ 


( 2 ) 


where fi,v = x,y, z are the components of the electric 
field, Kc is a proportionality parameter that is deter¬ 
mined numerically, Vp is the speed of the emergent 17(1) 
gauge photon, and /? is the inverse temperature. In the 
zero temperature limit —>■ oo), the above correlation 
turns into the 1/r^ dipolar-like correlation in real space 
|4] (see Fig. [^). 

In Fig. 1^ we depict our numerical results for the elec¬ 
tric field correlation {Ez{k)Ez{—k)) along different mo¬ 
mentum cuts of Brillouin zone. We have chosen the path 
k = (|, g, |) in Fig. to fit our data against Eq. [^and 
obtain the photon vdocity Vp = 0.0020(5)a, where a is 
the FCC lattice constant. The curve in Fig. Ep is plot of 
Eq. [fusing the fitted parameters from Fig.|^, along the 
path k = (0,q, 0). The fits and data match prefectly in 
Fig. and b and it provides a unbiased test of the data 
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FIG. 4. (Color online). Electric field correlators at a direc¬ 
tion. Dashed curves are fitted theoretical results using Eq. 
We choose tjV = 0.055, p/E = 1 inside the CMI-1/4 with 
system size L = 12 and inverse temperature /3 = 1000. 


as only one set of parameters is used to fit the data along 
different momentum cuts. 

Due to the gapless gauge photon with dispersion Ck = 
Uplkj in the Coulomb liquid phase, we expect that the 
heat capacity C'(T) « T^7r^a^/(15t;p) for the CMI-1/4 
at zero temperature limit, where a is the FCC lat¬ 
tice constant and is the photon speed. As shown 
in Fig. [^, this dependence of the low-temperature 
heat capacity is clearly obtained for the CMI-1/4 for 
T < 0.002E, and we fit to the theoretical prediction and 
obtain Vp = 0.0017(4)a. As there is no continuous sym¬ 
metry breaking and thus no gapless Goldstone mode in 
the CMI-1/4, the heat capacity is a strong indica¬ 
tion of the presence of linearly dispersive gapless gauge 
photon m- Moreover, the photon speed extracted from 
heat capacity is consistent with the one obtained previ¬ 
ously from electric field correlation. These can be served 
as two independent evidences that the CMI-1/4 is indeed 
well described by the compact 17(1) QED. In Fig. [^, we 
depict the finite temperature behavior of the entropy for 
the CMI-1/4 phase |20j . 

Discussion. All exisiting theoretical and numerical 
works except Ref. [T5I consider the CMI-1/2 state of the 
hard-core boson model (or quantum spin ice in the spin 
language) on the pyrochlore lattice. We have carried out 
the first numerical study of the full phase diagram of 
the model and provide the numerical evidence for the 
bosonic CMIs in the phase diagram. Our numerical re¬ 
sults of Coulomb liquid phases of the bosonic CMIs are 
relevant to the charge sector physics of electron CMIs in 
the pyrochlore materials mm- One future direction is to 
include disorders and to understand if any intermediate 
Bose-glass-like phase occurs between the CMIs and the 
superfluid. If such a phase do exist for a disordered sys¬ 
tem, the physical property might be very different from 
the Bose glass for the conventional boson Hubbard model 


FIG. 5. (Golor online). Temperature dependence of (a) the 
specific heat C and (b) the entropy S for t/V — 0.055, fj./V = 
1 inside the GMI-1/4. The (magenta) dashed lines in (a) is 
the fit to specific heat of Coulomb liquid phase. The (orange) 
dotted-dash lines highlight the crossover temperature Do US]- 
The (green) horizontal lines denotes the Panling entropy of 
classical spin ice Sp = R| ln(|), and the entropy of free spins 
S = Rln(2) at the high temperature limit. 


with disorders m- 
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